{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "897bb39b-547f-4993-a49c-2018cc9648e3",
   "metadata": {},
   "source": [
    "# Backprojection image formation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a958b90-c742-4fa2-aef9-52202b7c8aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torchbp\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import hamming"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d247ef25-faf9-469d-bf68-54c716f41a37",
   "metadata": {},
   "source": [
    "Use CUDA if it's available, CPU otherwise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4062655e-b0dd-4f95-ad06-493ffb13445a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if torch.cuda.is_available():\n",
    "    device = \"cuda\"\n",
    "else:\n",
    "    device = \"cpu\"\n",
    "print(\"Device:\", device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d43373-ef7f-4895-852b-d43941fe1201",
   "metadata": {},
   "source": [
    "Constant definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a31e5901-925c-49bb-839c-b1bbaabf9294",
   "metadata": {},
   "outputs": [],
   "source": [
    "nr = 128 # Range points\n",
    "ntheta = 128 # Azimuth points\n",
    "nsweeps = 128 # Number of measurements\n",
    "fc = 6e9 # RF center frequency\n",
    "bw = 100e6 # RF bandwidth\n",
    "tsweep = 100e-6 # Sweep length\n",
    "fs = 1e6 # Sampling frequency\n",
    "nsamples = int(fs * tsweep) # Time domain samples per sweep\n",
    "\n",
    "# Imaging grid definition. Azimuth angle \"theta\" is sine of radians. 0.2 = 11.5 degrees.\n",
    "grid_polar = {\"r\": (90, 110), \"theta\": (-0.2, 0.2), \"nr\": nr, \"ntheta\": ntheta}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bed5dbac-3f81-4a93-a47f-72b9b5fcaa39",
   "metadata": {},
   "source": [
    "Define target and radar positions. There is one point target at 100 m distance and zero azimuth angle.\n",
    "For polar image formation radar motion should be in direction of Y-axis.\n",
    "If this is not the case positions should be rotated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b18c7f-303b-44e1-b4ec-e56eaea7fff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_pos = torch.tensor([[100, 0, 0]], dtype=torch.float32, device=device)\n",
    "target_rcs = torch.tensor([[1]], dtype=torch.float32, device=device)\n",
    "pos = torch.zeros([nsweeps, 3], dtype=torch.float32, device=device)\n",
    "pos[:,1] = torch.linspace(-nsweeps/2, nsweeps/2, nsweeps) * 0.25 * 3e8 / fc\n",
    "\n",
    "# Not used in this case\n",
    "vel = torch.zeros((nsweeps, 3), dtype=torch.float32, device=device)\n",
    "att = torch.zeros((nsweeps, 3), dtype=torch.float32, device=device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4047b6f-d590-4ee9-a4fc-c3e13f7d770c",
   "metadata": {},
   "source": [
    "Generate synthetic radar data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73133328-4cdb-4cbb-b4cd-ac2ce81301ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Oversampling input data decreases interpolation errors\n",
    "oversample = 3\n",
    "\n",
    "with torch.no_grad():\n",
    "    data = torchbp.util.generate_fmcw_data(target_pos, target_rcs, pos, fc, bw, tsweep, fs)\n",
    "    # Apply windowing function in range direction\n",
    "    w = torch.tensor(hamming(data.shape[-1])[None,:], dtype=torch.float32, device=device)\n",
    "    data = torch.fft.fft(data * w, dim=-1, n=nsamples * oversample)\n",
    "\n",
    "data_db = 20*torch.log10(torch.abs(data)).detach()\n",
    "m = torch.max(data_db)\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(data_db.cpu().numpy(), origin=\"lower\", vmin=m-30, vmax=m, aspect=\"auto\")\n",
    "plt.xlabel(\"Range samples\")\n",
    "plt.ylabel(\"Azimuth samples\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963bda9b-1f72-48f6-9ee1-d0ea54ff8cb9",
   "metadata": {},
   "source": [
    "Image formation.\n",
    "Hamming window was applied in range direction so low sidelobes in range are expected.\n",
    "Azimuth direction has no windowing function and high sidelobes (Highest -13 dB) are expected.\n",
    "Azimuth sidelobes could be decreased by windowing the input data also in the other dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd6b005-5d5f-432f-ba40-b4db21699609",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_res = 3e8 / (2 * bw * oversample) # Range bin size in input data\n",
    "\n",
    "img = torchbp.ops.backprojection_polar_2d(data, grid_polar, fc, r_res, pos, vel, att)\n",
    "img = img.squeeze() # Removes singular batch dimension\n",
    "\n",
    "img_db = 20*torch.log10(torch.abs(img)).detach()\n",
    "\n",
    "m = torch.max(img_db)\n",
    "\n",
    "extent = [*grid_polar[\"r\"], *grid_polar[\"theta\"]]\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(img_db.cpu().numpy().T, origin=\"lower\", vmin=m-30, vmax=m, extent=extent, aspect=\"auto\")\n",
    "plt.xlabel(\"Range (m)\")\n",
    "plt.ylabel(\"Angle (sin radians)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3d87a9-7c0b-4f46-bf05-873db9d7a336",
   "metadata": {},
   "source": [
    "Image entropy. Can be used as a loss function for optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b44fa95-3a4e-4b95-900b-3a89a8290d06",
   "metadata": {},
   "outputs": [],
   "source": [
    "entropy = torchbp.util.entropy(img)\n",
    "print(\"Entropy:\", entropy.item())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "768dc700-737b-4787-924d-a734ff9e7832",
   "metadata": {},
   "source": [
    "Convert image to cartesian coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d5ab162-2b7f-43d9-af03-1e23eb55f8db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Origin of the polar coordinates\n",
    "origin = torch.zeros(2, dtype=torch.float32, device=device)\n",
    "# Cartesian grid definition\n",
    "grid_cart = {\"x\": (90, 110), \"y\": (-10, 10), \"nx\": 128, \"ny\": 128}\n",
    "\n",
    "# Complex image can also be interpolated, but takings absolute value first is faster\n",
    "img_cart = torchbp.ops.polar_to_cart_linear(torch.abs(img), origin, grid_polar, grid_cart, fc, rotation=0)\n",
    "\n",
    "img_db = 20*torch.log10(img_cart).detach()\n",
    "\n",
    "m = torch.max(img_db)\n",
    "\n",
    "extent = [*grid_cart[\"x\"], *grid_cart[\"y\"]]\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(img_db.cpu().numpy().T, origin=\"lower\", vmin=m-30, vmax=m, extent=extent, aspect=\"equal\")\n",
    "plt.xlabel(\"Range (m)\")\n",
    "plt.ylabel(\"Cross-range (m)\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
